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Abstract. A two-step shape reconstruction method for diffuse optical to- 
mography (DOT) is presented which uses adjoint fields and level sets. The 
propagation of near-infrared photons in tissue is modeled by the time- 
dependent linear transport equation, of which the absorption parameter 
has to be reconstructed from boundary measurements. In the shape re- 
construction approach, it is assumed that the inhomogeneous background 
absorption parameter and the values inside the obstacles (which typically 
have a high contrast to the background) are known, but that the num- 
ber, sizes, shapes, and locations of these obstacles have to be reconstructed 
from the data. An additional difficulty arises due to the presence of so- 
called clear regions in the medium. The first step of the reconstruction 
scheme is a transport-backtransport (TBT) method which provides us with 
a low-contrast approximation to the sought objects. The second step uses 
this result as an initial guess for solving the shape reconstruction problem. 
A key point in this second step is the fusion of the 'level set technique' for 
representing the shapes of the reconstructed obstacles, and an 'adjoint-field 
technique' for solving the nonlinear inverse problem. Numerical experi- 
ments are presented which show that this novel method is able to recover 
one or more objects very fast and with good accuracy. 



1. Introduction 



1.1. A novel two-step shape reconstruction scheme. In a recent paper |12, 13| 
a two-step shape reconstruction method was introduced for the situation of 2D 
cross-borehole electromagnetic (EM) imaging. In that paper, electromagnetic cross- 
borehole tomography was modeled as an inverse problem for the 2D Helmholtz 
equation. By reformulating the inverse problem as a shape reconstruction problem, 
and applying a novel level set based shape reconstruction scheme, it was possible to 
recover nontrivial shapes from noisy, limited-view, multifrequency data. 



The method presented in |12, 13| is quite general. In the present paper, we 



demonstrate how this novel two-step shape reconstruction scheme can be applied 
to a completely different imaging situation. We are now interested in the imaging 
of absorption characteristics in human tissue from near-infrared light data. This 
imaging technique has become known as Diffuse Optical Tomography (DOT), see 
for example In contrast to the majority of approaches in DOT which use a 
diffusion approximation to the linear transport equation for the inversion task, we 
consider here the (time-dependent) linear transport equation itself, which is generally 
accepted as being the correct model for the photon propagation in tissue. 



This article was written during the stay of the author at the Mathematical Sciences Research 
Institute (MSRI), Berkeley, from October to December 2001. The research presented here is sup- 
ported in part by the NSERC CRD Grant 80357, and research during the stay at MSRI is supported 
by NSF grant DMS-9810361. 



2 



OLIVER DORN 



The inversion method is a two-step shape reconstruction scheme which tries to 
recover the unknown shapes of abnormahties in the absorption coefficient in the 
tissue. The first step consists of an inverse scattering algorithm for DOT which is 
called a 'Transport-Backtransport' (TBT) method and which was originally intro- 
duced in |l^. In the situation of DOT, this TBT algorithm yields already during 
the early iterations a relatively good approximation to the locations of the unknown 
obstacles, but requires many more iterations in order to get better estimates for the 
correct extension and shapes of these objects. Therefore, the result of these early 
sweeps is only used as an initial guess for initiating the second step of the two-step 
scheme, which is a level set based shape reconstruction method. This second step 
calculates successive corrections for the shapes of the obstacles which aim to reduce 
the misfit in the data. Here, the level set representation of the shapes is chosen 
because it enables us to easily describe and keep track of complicated geometries 
which might arise during the evolution process. 

Both, the initialization step as well as the level set based shape evolution step 
employ an 'adjoint-field technique' for the inversion which has the very useful prop- 
erty that the inverse problem can be solved approximately by making two uses of 
the same forward modeling code. Using a somewhat oversimplified description of 
our technique, the updates to the level set function are obtained by first making 
one pass through the code using the parameter distribution (which is here the ab- 
sorption parameter) corresponding to the latest best guess of the level set function, 
and then another pass with the adjoint operator applied to the differences in com- 
puted and measured data. Then the results of these two calculations are combined 
to determine updates to the level set function. The resulting procedure is iterative, 
and can be applied successively to parts of the data, e.g., data associated with one 
source location can be used to update the model before other source locations are 
considered. Similar techniques have been successfully applied recently as part of 
iterative nonlinear inversion schemes in ||, 10 , [ill, 32, 44|. For alternative inversion 
schemes in DOT see for example fl], [l^, Oj , gO], |25|, ^ and the references 

therein. 



1.2. The level set method. The level set method was originally developed by 
Osher and Sethian for describing the motion of curves and surfaces |36, 41 1. Since 
then, it has found applications in a variety of quite different situations. Examples are 
image enhancement, computer vision, interface problems, crystal growth, or etching 
and deposition in the microchip fabrication. For an overview see [p2|| . 

The idea of using a level set representation as part of a solution scheme for in- 
verse problems involving obstacles was first suggested by Santosa in ||4^. In that 
paper, two linear inverse problems, a deconvolution problem and the problem of re- 
constructing a diffraction screen, are solved by employing an optimization approach 
as well as a time evolution approach using the level set technique for describing the 
shapes. Santosa also outlines in that paper how the two presented methods can be 
generalized to nonlinear shape recovery problems. More recently, a related method 



was applied to a nonlinear shape recovery problem by Litman et al. in |25]. In 
that work, an inverse transmission problem in free space is solved by a controled 
evolution of a level set function. This evolution is governed by a Hamilton-Jacobi 
type equation, whose velocity function has to be determined properly in order to 
minimize a given cost functional. 
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In Dorn et al. [0, ^] a two-step shape reconstruction algorithm was presented 
for the situation of electromagnetic cross-borehole tomography in 2D. In contrary 
to the work mentioned above, this approach is not based on a Hamilton-Jacobi 
type equation, but uses an optimization approach instead. Moreover, a very specific 
inversion routine, an adjoint-field technique, is employed which calculates an update 
for the most recent level set function directly by just solving one forward and one 
adjoint problem on the computer. The present paper is based on this approach, and 
a detailed discussion will be given in the following sections. 

As further examples for related and very interesting recent approaches we also 
want to mention here the work described in Burger Elangovan et al [H], Hin- 



termiiller et al |21|, Ito et al [22], Jackowska-Strumillo et al |23[|, Osher et al [37|, 



Ramananjaona et al. [|38[, Sethian et al [43|, and Zhao et al 



The paper is organized as follows. In section g we will present the basic equations 
of DOT in a form convenient for development of the shape reconstruction technique. 
In section ^ we formulate the shape reconstruction problem and introduce the level 
set formulation of this problem. In section ^, we derive the basic shape reconstruction 
algorithm using level sets and adjoint fields. Section |5| describes shortly the proce- 
dure employed here for finding a good starting guess for our shape recontruction 
algorithm. In Section ^ numerical experiments are presented which demonstrate the 
performance of the algorithm in different situations. The final section summarizes 
the results of this paper and indicates some directions for future research. 

2. The physical experiment in DOT 

2.1. The linear transport equation. We model the propagation of photons in 
the medium by the time-dependent linear transport equation 

+ e-\/u{x,e,t) + {a{x) +b{x))u{x,e,t) - b{x) / r]{e ■9')u{x,e' ,t)d9' 
ot Js"-'^ 

= q{x,9,t) in Q x S''~^ x [0,T] (1) 

with initial condition 

u{x,9,0)=0 in 17x5""-^ (2) 

and boundary condition 

u{x,9,t) = on r_. (3) 

Here, 

r± := {{x,9,t) GdnxS'""-^ X [0,r], ±iy{x)-9 > o}. 

O is a convex, compact domain in R", n = 2,3, with smooth boundary dO,. In 
our numerical experiments, we will only consider the case n = 2, but the algorithm 
extends in a straightforward way to n = 3. v^x) denotes the outward unit normal 
to 50 at the point x £ dQ, and u{x, 9, t) describes the density of particles (photons) 
which travel in Q at time t through the point x in the direction 9. The velocity c of 
the particles is assumed to be normalized to c = lcms~^ and has been neglected in 
the formulation of (|l|). 

a{x) is the absorption cross-section (in short 'absorption'), b{x) is the scattering 
cross-section, and /u(x) := a{x) + b{x) is the total cross-section or attenuation. These 
parameters are assumed to be real, strictly positive functions of the position x. The 
quantity is the mean free path of the photons. Typical values in DOT are 



4 



OLIVER DORN 



a » 0.1 - 1.0cm~\ b ^ 100.0 - 200.0cm"\ /i"^ w 0.005 - 0.01cm, see for example 
|24|. The scattering function r]{6 -O ) describes the probabihty for a particle entering 



a scattering process with the direction of propagation 6 to leave this process with 
the direction 0. It is normalized to 

r]{e-e')de = i, (4) 



5" 

which expresses particle conservation in pure scattering events. The dot-product in 
the argument indicates that 77 depends only on the cosine of the scattering angle 
COST? = 6-6 and, in particular, is independent of the location of the scattering 
event x. These latter assumptions simplify the following calculations, but the recon- 
struction method can be extended to more general scattering functions as well. In 
our numerical experiments, we will use the following Henyey-Greenstein scattering 



function [24| 



2(1 + g'^ — Ig cos ■uy^i^ 

with — 1 < < 1. The parameter g in (^) is the mean cosine of the scattering 
function. Values of g close to one indicate that the scattering is primarily forward 
directed, whereas values close to zero indicate that scattering is almost isotropic. In 
our numerical experiments, we will choose g to be 0.9, which is a typical value for 
DOT. 

The initial condition (Q) indicates that there are no photons moving inside of 
at the starting time of our experiment. The boundary condition (|3|) indicates 
that during the experiment no photons enter the domain i7 from the outside. All 
photons inside of O originate from the source q which, however, can be situated at 
the boundary 50. We mention that an alternative way of describing laser sources 
at the boundary is to use an inhomogeneous boundary condition instead of (^) and 
a zero source (7 = 0. Both choices are equivalent, and the derivation of the TBT- 
method as well as the shape reconstruction method are almost identical for these 
two possibilities. 

We consider the problem d!])-® for p different sources qj, j = 1, . . . ,p, positioned 
at do,. Typical sources in applications are delta- like pulses transmitted at time 
t = at the position sj € di} into the direction Oj, which can be described by the 
distributional expressions 

qj{x,e,t) = 6,{sj)6t{0)Se{9j), j = 1, . . . ,p, (6) 

with i^{sj) ■ 9j < 0. We assume that a given source qj gives rise to the physical fields 
Uj{x,6,t) which are solutions of Our measurements consist of the outgoing 

flux across the boundary Oil which has with the form 

Gj{x,t) = f u{x) ■Ouj{x,e,t)de on dQ x [0,T], (7) 

for j = 1, . . . ,p. With these definitions, we are now able to formulate the inverse 
problem which we will consider in the following. 

Definition 2.1. Inverse problem of DOT: Assume that we measure for p different 
sources qj, j = 1, . . . ,p, the data Provided with this information and knowing 
the scattering functions ij and b, how can we reconstruct the coefficient a{x) of ^ 
inside of Vt? 
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In general, also the scattering parameter b{x) is part of the inversion problem. 
However, for simplicity, we will restrict ourselves in this paper to the determination 
of only one parameter, namely a{x), and will assume that the scattering parameter 
b{x) is given. Extensions of this method will be addressed in our future research. 



2.2. Clear regions. An additional difficulty arises in our application of DOT due 
to the presence of so-called 'clear regions' in the imaging domain. These are regions 
where the scattering coefficient b{x) is very small or zero, such that the dominant 
behavior of the traveling photons inside these regions is transport rather than dif- 
fusion. Most of the imaging methods in DOT are purely based on the diffusion 
approximation to the linear transport equation, such that they cannot be used in 
applications where clear regions are present. Clear regions play an important role 
for example when imaging the head of neonates or in functional brain imaging. See 
|, 

Although we assume in this paper that the locations of the clear regions are known 
a priori (which in realistic applications is not always the case such that the shapes 
of the clear regions become part of the imaging task themselves), the presence of 
these regions complicates the numerical inversion process. 

One problem caused by the clear layers is their effect on the sensitivity structure 
of the data with respect to the unknown parameters. The clear regions (and, in 
particular, the clear layers which are surrounding the region of interest, see figure 
I) act like 'wave guides' for the photons, since the propagation is much faster in- 
side these layers than it is in the diffusive regions. Therefore, the photons which 
have travelled through these clear layers dominate in particular the early parts of 
the time-dependent measurements. Unfortunately, these photons do not carry much 
information related to the internal structure of the medium. In absence of clear lay- 
ers, the early photons are usually assumed to yield the reconstructions with the best 
resolution. If clear layers are present, however, these early photons are now mixed 
with those photons that travelled through the clear layers, such that it is difficult 
to extract the necessary information about the interior domain from the early data. 
Therefore, we will in our numerical reconstructions only use photons which arrived 
at later times at the receiver locations. Our sensitivity analysis indicates that these 
later photons are less affected by the clear layers. We will show that the informa- 
tion of the later arriving photons is sufficient for achieving good reconstructions of 
the absorption distribution inside the regions of interest, which might however have 
a slightly reduced resolution. See also the more detailed discussion of sensitivity 
functions in section 4.7. 



3. The shape reconstruction problem 

In this section we formulate the shape reconstruction problem which we want to 
solve, and cast it in a form which makes use of the level set representation of the 
domains. 



3.1. The shape reconstruction problem in DOT. To start with we introduce 
some terminology which we will use throughout the paper. 

Definition 3.1. Let us assume that we are given a constant a > 0, and a bounded 
function ai^ : ^ R. We call a pair {D,a), which consists of a compact domain 
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D CC and a bounded function a : 17 ^ R, admissible if we have 

a|z) = a, «b\D = abb\L>- (8) 

In other words, a pair {D, a) is admissible if a is equal to a preassigned constant 
value a inside of D, and equal to the preassigned background absorption Cb outside 
of D. The domain D is called the scattering domain. 

Remark 3.1. For an admissible pair {D,a), and for given a, a^, the absorption a 
is uniquely determined by D. 

With this definition, we can now formulate the shape reconstruction problem. 

Shape reconstruction problem. Let us assume that we are given a constant 
a > 0, a bounded function Cb : 17 — > R, and some data Gj as in Find a domain 
D such that the admissible pair {D, a) reproduces the data, i.e. (|^) holds with uj 
being a solution of ([l|)-@ with a = a for j = 1, . . . ,p. 

Using the same notation and assumptions as in definition 3.1, we want to formu- 
late another inverse problem which we will call the inverse scattering problem and 
which will play an important part when solving the shape reconstruction problem. 

Inverse Scattering Problem. Let us assume that we are given a bounded func- 
tion ttb : O ^ R, and some data Gj as in (|^. Find a bounded function : ^ R 
such that a = Cb + Og reproduces the data, i.e. holds with Uj given by (|l|)-(§) 
and with a = d for j = 1, . . . ,p. 

The inverse scattering problem gives rise to the decomposition 

o = Cb + in (9) 

In other words, the absorption distribution a is decomposed into the background 
distribution Ob and the perturbation which we will refer to in the following as the 
obstacle. 

Solving the shape reconstruction problem requires only to find the shape of the 
domain D, since the function a is then uniquely determined by (P). Solving the 
inverse scattering problem, on the other hand, amounts to finding the entire function 

from the given data, which is much harder to do. However, it will turn out that 
finding a good approximate solution of the inverse scattering problem is much easier 
and faster to achieve and will provide us with a very good initial guess for starting 
our shape reconstruction routine. 



Definition 3.1 allows us to formulate a first version of the strategy which we want 



to use for solving the shape reconstruction problem. 



Strategy for solving the shape reconstruction problem. Construct a series of 
admissible pairs {D^"\ o^"^); n, = 0, 1, 2, . . . , such that the misfit between the physi- 
cally measured data and the calculated data corresponding to {D^^\ a^"^) decreases 
with increasing n, and ideally, i.e. in absence of noise, tends to zero in the limit 
n — > oo. Use the approximate solution of the inverse scattering problem to initialize 
this series by determining a good starting element {D'^^\a^^^). 
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3.2. The domains D^^K In our numerical examples, each of the domains Z)*-"-* 
which we are looking for can be given as a collection of a finite number L„ of 

(n) 

disjoint, compact subdomains ,1 = 1,... , Ln, with 

jjin) ^ y ^(n)^ jjin) ^ ^(n) ^ g ^ ^^^^ 

The shapes of these subdomains -D;"^ can in principle be arbitrary. In particular, 
they are allowed to be multiply connected, and to enclose some 'cavities' or 'holes' 
filled with background material, and moreover the number Ln of these subdomains 
might (and usually does) vary with the iteration number n. For the derivation of 

(n) 

the inversion method, we assume that the boundaries dD^ of these domains are 
sufficiently smooth. 

It is essential for the success and the efficiency of the reconstruction scheme to 
have a good and flexible way of keeping track of the shape evolution during the 
reconstruction process. The method we have chosen in our reconstruction algorithm 
is a level set representation of the shapes as it was suggested by Santosa [|^. This 
representation has the advantage that the level set functions, which are in principle 
only used for representing the shapes, can in a natural way be made part of the 
reconstruction scheme itself. Doing so, it is not necessary anymore to refer to the 
shapes of the domains until the reconstruction process is completed. The final shape 
is then recovered from the representing level set function easily. In the following we 
will discuss in a more formal way how this can be achieved. 

3.3. Level set representation of the domains L'^"^ Assume that we are given 
a domain D CC 0. The characteristic function xd '■ ^ {0, 1} is defined in the 
usual way as 

= [l ; lln\D. ^^1) 

Definition 3.2. We call a function (p : 17 ^ R a level set representation of D if 

Xd{x) = -^4>{x) on Q. (12) 
where : $7 — > {0, 1} is defined as 



For each function : ^ R there is a domain D associated with (p by (|T^,(13) 
which we call D[(j)]. It is clear that different functions (l)i,4>2, 4>i 7^ (1)2, can be as- 
sociated with the same domain -D[(/>i] = D[(p2], but that different domains cannot 
have the same level set representation. Therefore, we can use the level set represen- 
tation for unambiguously specifying a domain D by any one of its associated level 
set functions. 

The boundary F = dD[(j)\ of a domain -D[0], represented by the level set function 
(j), is defined as 

F = {x E 11 : for all p > we can find xi,X2 € Bp{x) (14) 

with (/)(xi) > and (/)(x2) < } 
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Definition 3.3. We call a triple {D, a, (p), which consists of a domain D CC and 
bounded functions a,(f) : $7 — s- R, admissible if the pair (D, a) is admissible in the 
sense of definition 3.1, and (j) is a valid level set representation of D. 

Remark 3.2. For an admissible triple {D, a, (p), and for given d, a^, the pair (D, a) 
is uniquely determined by (j). 

We use these definitions to reformulate our shape reconstruction problem. 

Level set formulation of the shape reconstruction problem. Given a constant 
a > 0, a background distribution Ob, and some data G* as in Find a level set 
function (j) such that the corresponding admissible triple {D,a,(l)) reproduces the 
data, i.e. (^) holds with Uj being a solution of (|l])-@ with a = a for j = 1, . . . ,p. 

The strategy for solving this shape reconstruction problem has to be reformulated, 
too. It reads now as follows. 



Strategy for solving the reformulated shape reconstruction problem. Construct 
a series of admissible triples {D^^\ a^^\ </'^"^); n = 0, 1, 2, ... , such that the misfit 
between the data (|^) and the calculated data corresponding to (D^"), a^"^ 
decreases with increasing n, and ideally, i.e. in absence of noise, tends to zero in 
the limit n — > oo. For finding this series we only have to keep track of a*-"^ and 
(h^'^\ but not of D^^\ The function a(") is needed in each step for solving a forward 
problem (|l|)-@, and a corresponding adjoint problem. The knowledge of 0^"^ is 
used in each step to determine a*-"-*. The final level set function (f)'^^\ which satisfies 
some stopping criterion, is used to recover the final shape Z)*^^) via ([l^. 

4. Solving the shape reconstruction problem 

In this section we derive the basic shape reconstruction method which uses adjoint 
fields and the level set representation introduced above. The initializing procedure 
for this reconstruction routine will be discussed in section |5[ 

4.1. Function spaces. We want to specify now the function spaces which we will be 
working with. The main objective of this section is to introduce the inner products 
on these function spaces, which will become important when defining the adjoint 
linearized operators in the following sections. 

The natural spaces in the framework of the linear transport equation are L}- 
spaces for the fields u, and L°°-spaces for the parameters a and 6, due to the physical 
significance of the corresponding norms. This is because the L^-norm of a function 
u{x, 9, to) describes physically the number of particles which are in the system at 
a given time tQ. However, in our derivation of the shape reconstruction method it 
will be more convenient to work with L^-spaces instead, since deriving the adjoint 
operators used in the inversion routine is more straightforward when using these 
spaces. See for example |5|, ^, ^, |l^ for more information concerning L^-spaces in 
linear transport theory. The spaces which we will consider are 

Y = L'^inx S""-^ x[0,T]), U = L\nx S""-^ x[0,T]), (15) 



Z = L^idQ. X [0,r]), F = L^{n), $ = L^{n), (16) 
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equipped with the standard inner products. Here, Y is the space of sources q, U is 
the space of states u, Z is the space of measurements or data G, F is the space of 
parameters a, and $ is the space of level set functions (p. We mention that in fact 
our level set functions </» will be assumed to have a certain regularity, since in the 
derivation of the reconstruction scheme we will make use of first order derivatives of 
these functions. 

4.2. Operators. In the following, we will introduce some operators which will en- 
able us to formulate the shape reconstruction problem in a way suitable for deriving 
the inversion algorithm. 

Given a constant a and a bounded function a^, : ^ R. Then, with each level 
set function G <I> a uniquely determined 'absorption perturbation' or 'obstacle' 
A(0) is associated by putting 

With (jl^) we can write this also as 

A{^){x) = ^^{x){a-a^{x)) , x G fi. (18) 

Notice that the operator A is chosen such that the triple {D, a, (p) with a = ab + A(0) 
and domain D[(f)] forms an admissible triple {D, a, (f)) in the sense of definition |3.3| . 
Moreover, for {D, a, cp) an admissible triple, we see that is just the absorption 
perturbation 

A(0)(x) = a,(x) = XD{x){a - ai,{x)), xeQ. (19) 

Let us assume now that we are given a background absorption Ob and that we have 
collected some data Gj which correspond to the 'true' absorption distribution 

a = Ob + a^, (20) 

where is the 'true' obstacle. We consider the following measurement operators 

Gj : F — > Z, Gj{a,){x,t) := / u{x) ■ 6 Uj{x,e,t) dO, (21) 

j = 1,... ,p, where Uj solves ([l|)-(§) with the source qj G Y and the parameters 
a = Ob + ffls- For the 'true' parameters we ask (21) to coincide with the data Gj 

Gj{d,){x,t) = Gj{x,t) forj = l,...,p (22) 



on dil. X [0,T]. Assuming Gj G Z we want to determine such that (22) is valid. 
For given data Gj, j = 1, . . . ,p, we define furthermore the residual operators 

Rj : F — > Z, Rj{a,){x,t) = Gj{a,){x,t) - Gj{x,t). (23) 

From ( p2[ ) it follows that for the 'true' obstacle the residuals vanish, 

Rjid^) = forj = l,... ,p, (24) 

if the data are noise-free. 

The forward operators Tj which map a given level set function G $ into the 
corresponding mismatch in the data are defined by 

Tj-.^^ Zj , TjicP) = Rj{Am (25) 
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for J = 1, . . . ,p. The goal is to find a level set function G $ such that 

Tj{4>) = fori = l,...,p. (26) 
We mention that all three operators A, Rj and Tj are nonlinear. 

4.3. Linearized operators. For the derivation of the shape reconstruction algo- 
rithm, we will need expressions for the linearized operators corresponding to the 
nonlinear operators introduced above, and for their adjoints with respect to the 
given inner products. In this section, we define the linearized operators, and expres- 



sions for their adjoints are derived in section 4.5. 

Analogous to Santosa [^0| it is shown that, for a homogeneous background Ob, the 
infinitesimal response 6a{x) in the obstacle as{x) to an infinitesimal change (5(/)(x) of 
the level set function (j)(x) has the form 



da[x) = — [a — ObJ 



|V0(x)| 



(27) 

x€dD[<l>] 



The function 6a in (^) can be interpreted as a 'surface measure' on the boundary 
r = dD[(p]. Similar to (^), we define the linearized operator A'[(p] by 

{A'lm) (x) = -[a- a,{x)] ^^^^^rix) (28) 

where 5r{x) denotes the Dirac delta distribution concentrated on T = di}[(j)]. In 
this interpretation, (|27| ) describes an infinitesimal 'surface load' of absorption on T 
which has to be recovered from the mismatch in the data. 

However, the expression on the right hand side of (|28| ) is not an element of F 
which causes problems when we want to calculate the inner products defined in 



section iA. Therefore, we will replace the operator ( psj) by an approximation which 
maps from <I> into F and which will be more convenient for the derivation of the 
reconstruction method. 

For a given level set function S <I>, let F = dD[cj}] and Bp{T) = Uy^^Bpiy) a 
small finite width neighborhood of F with some given constant < /j <C 1. The 
(approximated) linearized operator A.' [(f)] is defined as 

m -.^^F, [mn) (^) = - [a - a,{x)] ^^Cp(r)XB,(r)(x) 

\^9\xj\ ^29) 



with XB,{v){x) given by (0). Here, Cp{T) = L{T) /Yo\{B p{T)) where L(F) = 
^Q^5Y{x)dx is the length of the boundary F, and Vol(i?p(F)) = JqXb {r){x)dx is 
the volume of Bp{T). For a very small p we will get a very large weight Cp(F), 
whereas for increasing p this weight Cp{T) decreases accordingly. The operator de- 
fined in ([29|) maps now from $ into F such that we can make use of the inner 



products defined on these spaces. 

We already mentioned that the term |V(/>(x)| in (27), (^), as well as the derivation 
of these expressions, implies some regularity constraint on (p. For example, (p G 
would be possible. Another possibility would be to use a 'signed distance function' 



as a standard representation of the boundary [42|. We do not want to specify the 
regularity of (j) at this point, but assume instead that it is 'sufficiently smooth' for 
our purposes. 
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The linearized residual operator R'j [a^] is defined by 

RAa,] : F — > Z, R.:[a,]5a = [ v ■ 6 vAx.O ,t) dO (30) 
where Vj solves the problem 

^ + e-vvj{x,e,t) + {a{x) + b{x))vj{x,e,t) - b{x) f 7]{9 ■e')vj{x,e' ,t)de' 

= -6a{x)uj{x,e,t) in J7 x 5""^ x [0, T] (31) 

with initial condition 

Vj{x,9,0)=0 in 0x5"-^ (32) 

and boundary condition 

Vj{x,d,t) = on r_. (33) 

Here, uj is the solution of (|l|)-(3) with source qj and with a = a^, + a^. Note that 
(|3l|) has the same structure as (1). The 'primary sources' g of (|l|) are now replaced 
by 'secondary sources' which are caused by the parameter perturbation 6a. The 
representation (^) can be derived by perturbing 

^ as + 6a , uj Uj + vj, (34) 

plugging this into (|^) and neglecting terms which are of higher than linear order in 
the perturbations 6a, Vj. A rigorous derivation of (^) in the framework of L^-spaces 
can be found in Q. 

As our third linearized operator, we introduce the linearized forward operator 
rj[0] by putting 

Tj[</.] -.^^Zj , T;[<p]6cj> = R'^ [A(,/.)] A'[cP]6cp. (35) 
All three operators A' [</>], i2^ [as], and are linear. 

4.4. A nonlinear Kaczmarz-type approach. The algorithm works in a 'single- 
step fashion' as follows. Instead of using the data (|^) for all sources simultaneously, 
we only use the data for one source at a time while updating the linearized residual 
operator after each determination of the corresponding incremental correction 6(/). 

To be more specific, let us assume that we are given a level set function (j)^"\x) 
and an obstacle a'f^\x) such that {D^'^^a^^ + a^"\ (/>'•")) forms an admissible triple 
in the sense of definition Using a data set Gj corresponding to the fixed source 
position qj, we want to find an update (5^^") to (f)^'^'^ such that for the admissible 
triple 

(z)('^+i),a, + a("+^),<A("+^)) := (36) 

+ + A(0(") + 5,/.(")), (/.(")+ 

the residuals in the data corresponding to this source vanish 

T,-(0("+i)) = T,-(0(") +(5</)(")) = 0. (37) 

Applying a Newton-type approach, we get from ( ^7] ) a correction for 0^") by 

solving 

rj[(/.(")]5</)(") = -T^""^) = - (UjUj-Gj) (38) 



12 



OLIVER DORN 



where uj satisfies (|l])-(|3|) with 

a(.) = aj.) + A(*(">)(.) = { ^j^j ; I I ^I*;"J1 ,„,j (39) 

Since we have only few data given for one source and one frequency, equation (^) 
usually will have many solutions (in the absence of noise) , such that we have to pick 
one according to some criterion. We choose to take that solution which minimizes 
the L^-norm of 50^") 

Min ||(50(")||2 subject to Tj(0("))(50(") = - (u^Uj - Gj) . (40) 

This solution can be formulated explicitly. It is 

= -t;[<^(")]* (t;[/")]t;[/")]*)"' {mju, - Gj) , (4i) 

where rj[(^'^")]* denotes the adjoint operator to rj[(/)(")]. 

After correcting </> by (/> — > ct) + 5(pj, where 54) j is given by (^), we use the updated 
residual equation (|3^) to compute the next correction 5(t)ji. Doing this for one equa- 
tion after the other, until each of the sources qj has been considered exactly once, 
will yield one complete sweep of the algorithm. This procedure is similar to the 
Kaczmarz method for solving linear systems, or the algebraic reconstruction tech- 
nique (ART) in x-ray tomography |pl| and the simultaneous iterative reconstruction 
technique (SIRT) as presented injR]. For related approaches in a variety of imaging 
problems we refer to |6|, [l^, 11, ^4 33, 34, 44]. 

4.5. The adjoint linearized operators. In order to calculate the minimal norm 
solution (41), we will need practically useful expressions for the adjoints of the 
linearized operators of section ( |4.3| ). We will present such expressions in this section. 
The calculation of the actions of these operators will typically require us to solve 
an adjoint transport problem. This explains the name 'adjoint-field method' of the 
inversion method employed here. 

To start with, a simple calculation gives us the following theorem. 

Theorem 4.1. The adjoint operator A' [(/>]* which corresponds to the linearized op- 
erator A'[0] is given by 

A'lcPY : F ^ $ , (A'[c/>]*<5a,) (x) = - [a - a,(x)] ^fi^C,(r)xB,(r)(^). 

The next theorem describes the adjoint operator i?j[as]* which corresponds to 
R'jla^]. Its proof can be found in |^], and is therefore omitted here. 

Theorem 4.2. Let z ^U* be a solution of the (adjoint) transport equation 
dz 



e ■Vz(x,e,t) + (a(x) +b(x))z(x,e,t) - b(x) T](9 ■e)z(x,e ,t)d9 

at 



gn 



= in nx 5"-^ X [0, T] (43) 

with 'initial' condition 

z{x,e,T)=0 on nxS""-^ (44) 
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and boundary condition 

z[r+ = z{x,t) := Cix,t). (45) 
Then the adjoint operator R'jla^]* is given by 

(r'M*C){x) = Ij{a,) := [ [ Uj{x,e,t)zj{x,e,t)dedt, (46) 

J[Q,T] JS"-i 

with Uj being a solution of with source qj and with a = Ob + Og. Notice that 

C{x,t) is applied uniformly into all directions 9 with u ■ 9 > Q. □ 



Finally, by combining theorems |4.l| and we get an expression for the adjoint 
operator Tj [</<]* which corresponds to the linearized forward operator Tj [</)]. It is 
described in the following theorem. 

Theorem 4.3. Let C, ^ Z . The adjoint operator Tj[(/)]* acts on C, in the following 
way 

T'^WC = A'[0]*ii;[A(c/>)]*C = i^^^^C,(r)xB,(r)(x)/,(A(0)) (47) 

where Ij{K{(j))) is given by and where Uj solves and Zj solves (^^-(^^ 

with replaced by A(0). 

4.6. Updating the level set function. Let us consider the operator Cj^n '■= 
Tj[0(")]Tj in (^) more closely. In order to calculate a correction 5</>mn by 



El|) we first have to apply the operator to the expression Mjuj — Gj. This 
operator is difficult to calculate, since typically a large number of forward problems 
has to be solved for this purpose if the background medium is inhomogeneous and 
the corresponding Green functions are not easily available. Moreover, the inversion 
of this operator is usually unstable due to the ill-posedness of the underlying inverse 
problem. Therefore, in our numerical experiments, we will usually just replace this 
operator by a suitably chosen approximate operator Cj which is easier to calcu- 
late and whose inversion is more stable. One possibility for such an approximation 
(which is adopted here) is to simply use a multiple of the identity, combined with an 
appropriate selection of receivers in the current step of the nonlinear Kaczmarz-type 
approach. We refer here for a more thorough discussion of this operator and its 
approximations to |ll|, 13, 34 1. With this approximation, we define 



:= C-'{m,u,-G,). (48) 
Next, in order to find Scp^^, we have to calculate Tj ]*Cj. Using we get 



(49) 



where Uj solves (|l|)-(^ with source qj and zj solves (p3|)-(^ with ( = Qj and 
a, = A(</.(")). 

In the following we want to introduce an approximation to the update formula 
( p9| ) which in our numerical experiments so far has shown to simplify and stabilize 
the reconstruction process. We mentioned already earlier that the level set represen- 
tation of a given contour is not uniquely defined. Therefore, it would be desirable 
to pick from all those possible level set functions representing a given contour that 
one with the most convenient properties. 
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Since we are updating the level set function in each step of our algorithm, it might 
happen that the denominator in ( ^9[ ) becomes very small, and numerical noise and 
roundoff-errors can destabilize the inversion routine. In order to avoid this problem, 
we rescale the level set function (/> after each update with the aim to keep |V(/)(x)| 
near the boundaries of the inclusions in a certain range of values. As a simple scaling 
rule, we multiply the level set function with a scaling factor such that the global 
minimum (or alternatively maximum) of the level set function is always equal to a 
predefined constant value. Notice that this rescaling procedure does not change the 
domains which are represented by the given level set function. In other words, we 
try to pick from the family of level set functions representing a given domain that 
one which has a convenient behavior of |V(/>(x)| in the neighborhood of the zero level 
set. 

We go one step further and make the assumption that, with this rescaling proce- 
dure, the gradient |V(/)(x)| along the boundary can be approximated in (|49| ) by some 
constant ci. Although, with this approximation, we are not marching anymore into 
the 'optimal' direction when calculating the updates, we will gain some stability in 
the inversion routine since we avoid dividing by small numerical derivatives (these 
are still possible even after the rescaling) in (|49| ) which are strongly influenced by 
unavoidable numerical noise and roundoff errors. Moreover, our numerical experi- 
ments show that these small deviations from the 'optimal' marching direction are 
easily corrected by the succeeding updates. 

With this modification, we get the following update formula for the level set 
function 

64>^^Hx) = - (a - a,(x))cr^ C,(r)xB,(r)(x)/,(A(0("))). (50) 

In ( ^0| ) we did not incorporate the rescaling procedure described above. 

Notice that, although we did not explicitly impose any regularity constraints 
on the updates (|50|), they are now in the range of i?^.[A(0("))]* (up to the factor 
a — ab(x)) which implicitly gives us some information about the regularity we can 
expect. Our numerical experiments so far indicate that the degree of regularity 
which is achieved by applying ( |50|) is typically sufficient 'for practical purposes' in 
our numerical examples. 

4.7. Sensitivity analysis. A valuable tool for designing an imaging experiment 
is a careful sensitivity analysis of the underlying inverse problem. The sensitivity 
functions which correspond to a nonlinear inverse problem are functions which de- 
compose the linearized residual operator (which is often called 'Frechet derivative') 



as well as its adjoint operator, as described in more details in |14, 15]. The sensi- 
tivity functions corresponding to one given source position and one given receiver 
location describe the sensitivity of this source-receiver pair to small variations of 
the parameters in the medium. Locations in the medium where the sensitivities are 
small have almost no influence on the measurements at this specific detector loca- 
tion. Therefore, hidden objects at these locations are difficult to recover unless other 
source-receiver pairs have larger sensitivity values at these specific locations. Ideally, 
we would like to design an experiment where certain source-receiver configurations 
have a sensitivity structure which is large on a very small area in the domain, and 
very small elsewhere. The corresponding data can then be used for determining 
the coefficients at these specific locations. If we find a set of sensitivity functions 
such that each location in the medium is 'probed' by such an ideal configuration. 
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we expect to get good and fast reconstructions. Unfortunately, in diffusion- type 
inverse problems like ours, the sensitivity functions are typically spread out over a 
relatively large area. In addition, they have a large dynamical range since they have 
very large values close to the source and receiver locations, and very small values in 
the regions of interest. Therefore, it is very difficult to probe localized areas inside 
the domain with suitable superpositions of sensitivity functions. The role of the 
sensitivity functions in our TBT-reconstruction scheme is described in the following 
theorem. 

Theorem 4.4. Transport sensitivity functions. Given a source We can 
find functions ^ai^ntr', Xsc) with the following property. 

1. ) 'Projection step': The action of the linearized residual operator R'la^] on the 
perturbation 5a in parameter space can be described as follows 

(^R'[a^]6a^ (Xrjtr) = j '^a{Xr,tr;Xsc)5a{Xsc)dXsc- (51) 

2. ) 'Backprojection step': The action of the adjoint linearized residual operator 
R [aj* on the vector ( in data space can be described as follows 

(R*C)iXsc)= / / '^a{Xr,tr;Xsc)CiXr,tr)dtrdXr. (52) 

^ ^ JdnJ[o,T] 

The functions ^a{xr,tr',Xsc) are called 'absorption transport sensitivity functions'. 
Here, Xr and tr denote receiver location and receiver time, respectively, and Xgc is 
an arbitrary scattering point in the medium. 

A derivation of this theorem, and a more detailed discussion of these sensitivity 
functions, can be found in [0]. We have displayed typical sensitivity functions 
for our imaging configuration in figure |^. The upper row of the figure shows the 
distribution of the scattering parameter for two different imaging situations. On the 
left, we have no clear layers, on the right we have a clear layer close to the boundary 
of the imaging domain which might represent for example the Cerebro-Spinal Fluid 
Layer (CFL) of the human head [||, 35, The source position considered here is 



indicated as point 'PI' in the figure on the left hand side. The receiver is located 
at the point 'P3'. (We are not considering here the also marked point 'P2'.) These 
points are the same in the situation shown on the right hand side. The second 
row displays the sensitivity functions with respect to the absorption coefficient for 
these two situations and for those data which are due to a source position at PI 
and taken at the receiver position P3 at an early time, namely i = 10 s (recall the 
normalization c = 1 cm/s). It can be seen that these early times are very sensitive 
to the clear layer. Our interpretation of this observation is that the early photons 
are very likely to have been channeled through this clear layer, and are therefore 
very sensitive to small variations of the absorption parameter inside the clear layer. 

The third row shows the same sensitivity functions, but at a later time t = 24 
s. We see that these later times are more sensitive to variations in the absorption 
parameter inside the domain and less to those inside the clear layer. Photons which 
arrive that late at the given receiver, are very likely to have travelled a long time 
in the scattering medium, and are therefore more affected by parameter variations 
inside the medium. However, we also observe that these later sensitivity functions 
are more 'spread-out' than those corresponding to earlier times, since the photons 
had more time to reach remote areas inside the imaging domain. 
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Table 1: The level set reconstruction algorithm. 

Initialization. 
n = 0; 

(Z)W,a(°),(/.W) given from TBT. 
Reconstruction loop. 

FOR i = 1 : Im perform Im sweeps 

FOR j = 1 : p march over source positions 

Cj = Cj^{MjUj - Gj); Uj solves (|l|)-(|) with a^, qj 

5<j)i'^) = -{a - ab(x))/j(A(0(")))xB,(r); Zj solves (H)-® with a(") and Cj 



^(n+i) ^ c^g) (</.(") + update level set function 

= Ob + A((/)("'+-'^)); n = n + 1; reinitialization n — > n + 1 

END 
END 

, Cb + A(</)*^"')), i;^**-"'-'); final reconstruction. 

Since the sensitivity structure of these later photons is more focused onto the do- 
main of interest, we will use only these data for solving our inversion task. The in- 
version process incorporates an iterative backprojection ('backtransport') step along 
these sensitivity functions, as indicated in theorem and explained in |14, 15 1. 



Since the functions along which the data are backprojected are very spread-out, we 
expect a loss in resolution compared to those reconstructions which can make use 
of the more focused earlier-time sensitivity functions. 

4.8. The algorithm. The nonlinear Kaczmarz-type method for shape reconstruc- 
tion using level sets is described in brief algorithmic form in table 1. Here, ?? is a 
relaxation parameter for the update of the level set function which is determined 
empirically. The operator CJ^ is taken to be simply a multiple of the identity, such 
that it can be considered as part of ry. The constant CpiV) could be calculated 
explicitly for the actual curve F^"^), or it could be approximated by some value cor- 
responding to a simple geometrical object (to give an example, in case of a single 
circle it would be CpiV) = {2p)~^). In our numerical experiments so far, however, it 
is simply considered as part of rj. The same holds true for ci. The constant p is in 
our numerical experiments chosen between one and two grid cells. The scaling factor 
Cls^ is determined after each update to keep the global minimum (or maximum) of 
the level set function at a constant value. 

Notice that we typically apply a very small relaxation parameter ry, which has 
the effect that the accumulated contributions of the earlier updates are not totally 
destroyed by the new correction, but instead have some 'inertia' which is the larger 
the smaller we choose ry. This slows down the whole reconstruction process to 
some degree, but, on the other hand, stabilizes it significantly. In our numerical 
experiments, we typically choose r/ such that the corrections to the boundaries of 
the shapes in each update are not much larger than one or two pixels in the normal 
direction to these boundaries. 
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Figure 1. Transport absorption sensitivity functions 

—^a{xr,tr',Xsc) for a source located at PI and a receiver at 
P3. Left column: homogeneous background. Right column: homo- 
geneous background with clear region close to boundary. Upper row: 
the domains. Middle row: at an early time step (tr = 10 s). Bottom 
row: at a later time step (tj. = 24 s). 

5. The TBT-method and the first guess 

For starting our shape reconstruction method using level sets we will need an 
initial guess {D^^\a^'^\(j)^'^^). Although it is possible just to create an arbitrary 
initial guess without using any data at all, there are several reasons for employing 
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100 120 



Figure 2. A typical data curve at one receiver in DOT. 



a preprocessing step which is relatively inexpensive and which yields a good initial 
guess for the level set reconstruction method. This preprocessing routine usually 
indicates already the approximate locations of possible inclusions, such that the 
amount of work which has to be done by the level set method is significantly reduced, 
and the algorithm converges much faster. 

We are employing an inverse scattering algorithm (the TBT-algorithm) for this 
purpose which was introduced in ||9, IC] and which shares basic features with the level 
set reconstruction scheme. Typically, in the diffusive regime, this inverse scattering 
algorithm yields already in the early sweeps good low-contrast approximations to the 
high-contrast objects. These low-contrast approximations usually indicate already 
the approximate locations of the unknown objects, but typically do not capture the 
correct number, sizes, shapes or contrasts of these objects. For this, many more 
iterations of the TBT scheme would be necessary. Therefore, it is advantageous 
to apply only a few steps of this inverse scattering method, and to switch then 
to a shape reconstruction formulation which explicitly takes into account the high- 
contrast assumption. We will show that in fact much better reconstructions can be 
achieved rapidly in these situations by the shape reconstruction method compared 
to the inverse scattering algorithm. We refer for a detailed description of the TBT- 
method to the publications ||9|, and describe in the following only the scheme 
which we use for extracting a suitable initial triple {D^^\a^'^\(l)^^^) from the TBT 
reconstructions aTBT(3;). 

Let as before D denote the sought domain. We assume that we are working in a 
high-contrast situation, such that exactly one of the following conditions is satisfied 



a 
a 



ab(x) » for all X £ D 
ab(x) <^ for all x £ D. 



(53) 
(54) 



Since we know a and ab(x), we know the constant 



(55) 
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Choose a threshold value < 7ls < 1 (in our numerical examples presented in 
section |6| we use 7ls = 0.9) and define 

Ols := 7LsmaxaTBT(a;). (56) 

For the level set zero Lq'^^ of (p^^^ we require that 

= {x e : aTBT(2;) = Cls} • (57) 

This means that we want all points of where the reconstruction a-Y-BTix) has exactly 
the value Cls to be mapped to zero by the level set function (j)'^^^ 

(j)^'^\x) = for all X G L^q\ (58) 
The level set function is now defined as 

(/>(0)(x) = C^°)sign(£») (a^s - a^B^{x)) , (59) 

where C^J is some suitably chosen scaling factor. Notice that (|58|) and (|5^) are 
consistent. The initial scattering domain D^^^ and the absorption a^^^ are defined 
as 

= D[^'^% = a, + A(0W). (60) 

Together with (f)^^^ they form an admissible triple {D^^\a^^\(j)^^^). 

6. Numerical Experiments 

6.1. The computational domain. In our numerical experiments we consider a 
domain of the size 5x5 cm^ which is divided into 50 x 50 pixels each having a size of 
1x1 mm'^. One time step lasts 0.2 seconds ([s]), and measurements are taken during 
the first 100 time steps or 20.0 s. Since the velocity of the light in the medium is 
normalized to c = 1.0 cm/s, each particle can travel 20 cm in that time. The 
directional variable is discretized into 12 direction vectors which are equidistantly 
distributed over the unit circle. The scattering law is a Henyey-Greenstein function 
(|5|) with g = 0.9. The space and time derivatives of (|l|) are discretized by using a 
finite difference scheme. A more detailed description of our numerical scheme for 
solving (ll)-(^ can be found in [P, [lO|. 

Our sources have a width of 5 pixels or 5 mm along the boundary. Using these 
slightly extended sources reduces artifacts which tend to appear close to the bound- 
aries, but is not essential for the presented method. In total, there are 16 of these 
sources, 4 at each side of the domain. These sources are emitting an ultrashort (one 
time step) laser pulse at time t = in the direction perpendicular to the boundary 
into the medium. They are arranged in a way such that only the central 20 mm of 
each side are illuminated during the whole imaging process. In other words, those 
parts of the boundary whose distance to any of the corners is less than 15 mm are 
not illuminated in our experiments. The reason for this source arrangement is the 
location of the clear layer close to the boundary, combined with the fact that our 
domain is rectangular rather than circular. By this source arrangement, we avoid 
that ballistic (i.e. unscattered) photons travel along the clear regions. 

Receivers are positioned at those parts of the boundary which have a distance 
(along the boundary) of at least 5 cm or 50 pixels to the source position. Therefore, 
only half of the boundary is actually covered by the receivers in a given experiment. 
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This arrangement has been chosen since the source-near receivers are not very sensi- 
tive to the regions of interest inside the domain. Therefore, they are not supposed to 
carry much information about the sought objects. In addition, we beheve that the 
approximation of the matrices Cj^n by a multiple of the identity is easier to justify 
with this source-receiver geometry, because of the reduced 'dynamic range' of the 
corresponding data. However, we do not have any proof for this assumption so far. 

Another important detail is that we are only using those photons as data which 
arrived at the receiver positions at times between t = 8 s and i = 20 s. In other 
words, the early photons are not used for the inversion task, since they are more likely 
to have been channeled through the clear region close to the boundary. We plan 
to use the information of these early photons in our future work for the additional 
determination of the structure of the clear layers from the measured data, which is 
here simply assumed to be known. 

The background of our imaging medium is homogeneous except of the embedded 
clear layers. These homogeneous regions have a scattering parameter of 6 = 100.0 
cm~^, and an absorption parameter of = 0.1 cm~^. Inside the clear regions, 
the parameters are b = 0.01 cm~^ and at = 0.01 cm~^. In the first two numerical 
experiments, we have only one clear layer close to the boundary and surrounding 
the interior imaging domain, see figures ||, ^, ^ In the third numerical experiment, 
we have two additional disc-shaped clear regions imbedded in the interior imaging 
domain, see figure ^. 



6.2. Three numerical experiments. The results of our first numerical experi- 
ment are shown in figure ^ We want to recover three disc-shaped objects in the 
interior of the domain from time-dependent data. The true configuration is shown 
in the bottom right image of the figure. The values inside the obstacles are a = 0.5 
cm~^ and b = 100.0 cm~^. The data are created by running our finite differences 
discretization of the linear transport equation on the correct model. So far, 

no noise is added to the data. We plan to investigate the effect of noise in the data 
more carefully in our future research. 

First, we run our preprocessing step (a TBT algorithm) in order to find a good 
estimate of the number and locations of the obstacles. The results of these runs 
after 5 sweeps and after 20 sweeps of the TBT method are shown in the first two 
images of the figure. Notice that one sweep of the TBT method consists of using 
the data for each of the sources exactly once, in a nonlinear Kaczmarz-type scheme. 
For more details of that scheme see section iA and We see that the 

locations and even the number of the unknown obstacles are already indicated in 
the reconstructions after 20 sweeps, but that the contrast and the shapes are far 
away from being correct. We use the result in order to start our level-set based 
shape reconstruction scheme. The figure shows the initial shape which is extracted 
from the TBT reconstruction, the reconstruction after six steps of the algorithm 
(which means after only using the information corresponding to six of the 16 source 
positions once), and after 16 steps or one sweep of the algorithm. We see that 
after only one additional sweep, now with the level set based shape reconstruction 
scheme, we have arrived at a very good estimate of the shapes and locations of the 
unknown obstacles. Additional iterations (sweeps) do not improve the reconstruction 
significantly, such that they are not displayed here. In this second step we have 
assumed that we know the correct contrast of the obstacles to the background. We 
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conclude that the exphcit incorporation of this additional information of the contrast 
in our reconstruction scheme enables us to get rapidly a much better estimate of 
the shapes of the obstacles than it would have been possible by just continuing the 
pixel-based TBT reconstruction scheme. 

Certainly, the final goal of the reconstruction algorithm will be to find the shapes 
as well as the contrast of the obstacles simultaneously from the given data. We 
will address this question in our future research. Notice that the shapes of the 
obstacles change topology during the evolution process. The level set method has 
no difficulties in keeping track of these topology changes. 

As mentioned, in practical applications we usually do not know the correct values 
of the absorption parameter inside the shapes. In our second numerical experiment, 
we investigate the behavior of the level set scheme in the case that we prescribe a 
slightly incorrect value for the contrast. Instead of = 0.5 cm^^ inside the obstacle 
we assume = 0.55 cm~^. The new configuration is displayed in the bottom right 
image of figure |^. The other images of the figure show the results of the algorithm. 
We see from these results that the effect of this model error on the reconstruction is 
very small. After only one additional sweep with the level set based reconstruction 
scheme, we arrive at a very good estimate for the shapes of the obstacles. These 
shapes also do not differ much from those where the correct value is used. That 
indicates that the data which are used here are not very sensitive to small changes 
in the parameter values. 

In our third numerical experiment we investigate a slightly more complex situ- 
ation. Now, there are two additional disc-shaped clear regions in the domain as 
shown in Figure |^. The bottom right image of the figure shows the true configu- 
ration. Although we assume that we know the boundaries of these additional clear 
regions, such that we do not have to recover them simultaneously with the unknown 
obstacles, these regions might have a negative impact on our ability to recover these 
obstacles due to the changed sensitivity structure. We refer to for visualizations 
of the sensitivity functions in a similar situation. As an additional difficulty, we 
assume now that each of the three objects which we are looking for has a different 
contrast to the background. The three values of the absorption parameter inside the 
objects are a = 0.4cm~^, a = 0.5cm~^, and a = 0.6cm~^. In our reconstruction 
process, we apply the constant value a = 0.4 cm~^ for all three unknown obstacles, 
such that two of the objects are reconstructed with an incorrect value. Figure |5| 
shows first the reconstruction of the TBT-method after 20 sweeps. This is used for 
extracting a first guess for our shape reconstruction scheme, as it is shown in the sec- 
ond image on the top row. The third image on the top row shows the reconstruction 
after 6 steps of the shape reconstruction scheme (which means that only the infor- 
mation corresponding to the first six source position has been used exactly once), 
after three sweeps, and after 10 sweeps of the scheme. As in our first two numerical 
examples, the algorithm delivered very good estimates of the shapes after only a 
few sweeps. The evolution of the norm of the residuals during the reconstruction is 
displayed in figure ^. The residuals converge now more slowly than in the first two 
examples, and the later sweeps still yield small improvements to the reconstructions. 
Notice also that the final reconstructions of the obstacles tend to be slightly more 
extended than the original objects in the case that a smaller contrast value was used 
for the inversion. This, of course, might have been expected since a larger area of 
absorption material is now necessary in order to have the same effect on the data 
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Figure 3. The correct value a = 0.5 cm~ is used for the recon- 
struction. Top row from left to right: Ctbt after 5 sweeps of TBT; 
after 20 sweeps of TBT; initial shape for level set algorithm. Bottom 
row from left to right: shape after 6 steps of level set algorithm; after 
16 steps (i.e. 1 sweep); true object. 
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Figure 4. An approximate value a = 0.55 cm~^ (instead of a = 
0.5 cm~^) is used for the reconstruction. Top row from left to right: 
Otbt after 5 sweeps of TBT; after 20 sweeps of TBT; initial shape 
for level set algorithm. Bottom row from left to right: shape after 6 
steps of level set algorithm; after 16 steps (i.e. 1 sweep); true object. 
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as the original object with a higher absorption value. Finally, we observe also here 
that the shapes do change topology in order to arrive at the final solution. This is 
handled easily by the level set formulation. 

7. Summary and future research directions 

We have presented a novel two-step shape reconstruction method for DOT which 
uses level sets and adjoint fields. The first step of the method employes the TBT 
method in order to determine a first guess for the shapes. Then, the second step 
uses this first guess to start a shape evolution scheme in order to complete the in- 
version task. The combination of these two schemes makes use of the fact that the 
TBT method typically finds a good guess for the locations of hidden obstacles in the 
early sweeps, and that the level set based shape reconstruction method is particu- 
larly fast once a good approximation for the shapes has been found. Therefore, the 
combination of these two methods enables us to get very fast and accurate recon- 
structions of the unknown objects. The reformulation of the inverse problem as a 
shape reconstruction problem can be considered as a regular ization technique. The 
incorporated prior information about the high contrast of the inclusions and the ap- 
proximately known values of the parameters stabilizes and speeds up the inversion 
task. We believe that also the simultaneous reconstruction of shapes and contrast 
will be possible with a slight modification of the algorithm, and plan to investigate 
this interesting problem in our future work. 

So far, we have been using mainly the information from later arriving photons 
for the inversion task, in order to cope with the specific structure of the sensitivity 
functions when clear layers are present. In our future work, the additional informa- 
tion of the early photons will be used in order to try to determine also the structure 
of the clear layers from the given data. Moreover, we plan to extend the imaging 
problem to finding the scattering parameter b{x) as well. An important question to 
be addressed is also the proper incorporation of additional regularization techniques, 
and the correct choice of the underlying function spaces. Especially the choice of the 
level set space is not obvious from the beginning. The proper choice of this space 
will ideally also be connected to our regularization strategy. We conclude that the 
novel two-step shape reconstruction technique presented here is extremely promising 
for the application in DOT, and that it poses many interesting and mathematically 
challenging questions which are open for future research. 
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Figure 6. Norm of the residuals during the level set reconstruction 
process shown in figure 
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